Distortion of the perfect lattice structure in bilayer graphene 
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We consider the instability of bilayer graphene with respect to a distorted configuration in the 
same spirit as the model introduced by Su, Schrieffer and Heeger. By computing the total energy 
of a distorted bilayer, we conclude that the ground state of the system favors a finite distortion. 
We explore how the equilibrium configuration changes with carrier density and an applied potential 
difference between the two layers. 

PACS numbers: 81.05.Uw 



I. INTRODUCTION 

A planar arrangement of carbon atoms covalently 
bound via sp 2 orbitals exhibits a honeycomb structure 
and is denoted graphene. Graphene rose rapidly to the 
forefront of research in condensed matter physics, mostly 
because of the peculiar electronic structure that emerges 
from its crystal structure, and the consequent wealth of 
rich and unexpected phenomena 1 . Seminal experiments 
on two dimensional crystals established graphene as an 
accessible reality 2 , and immediately unveiled numerous 
surprises, both on a fundamental level - like a new form 
of quantized Hall effect - and on a practical and tech- 
nological level - like the highly efficient field effect and 
high electronic mobility 3,4 . Most of the appealing phe- 
nomenology of graphene owes to the fact that electron 
dynamics in this system can be described in terms of 
chiral massless Dirac fermions 3 and, in fact, graphene 
does exhibit many properties characteristic of relativistic 
particles 5 ' 6 . 

Equally remarkable phenomena occur in bilayer 
graphene, which consists of two adjacent graphene planes 
stacked in the A-B fashion typical of graphite. Bi- 
layer graphene displays the same sample quality and 
quasi-ballistic transport characteristic of its single layer 
counterpart 7 , but brings also its share of new physics 
stemming from the nature of its charge carriers: chiral 
massive electrons 8 . Most interesting is the fact that, al- 
though gapless in its pristine form, a potential difference 
between the two layers opens a gap in the spectrum that 
can be controlled via chemical doping 9 or gating 10 ' 11,12 . 

Despite such favorable prospects, the amount of knowl- 
edge gathered in the context of bilayer graphene still lags 
behind the intensity committed to single layer graphene. 
In this brief paper, we address a particular aspect of the 
electron-phonon interaction in bilayer graphene, namely 



the tendency to relax the perfect crystal structure and 
generate a static, uniform deformation. This effect is in- 
spired, and similar in spirit, to the well known Peierls 
instability that occurs in polyacetylene chains. As shown 
by Su, Schrieffer and Heeger (SSH) 13 ' 14 , in polyacety- 
lene the one dimensional chain of carbon atoms has a 
half-filled electronic ground state that is unstable with 
respect to a spontaneous dimerization. This dimeriza- 
tion opens a sizeable gap in the spectrum that can be 
easily detected experimentally 15 . 

The instability we envisage for bilayer graphene is re- 
lated with the application of SSH's ansatz to the inter- 




FIG. 1: (color online) Transverse view (bottom) along the 
dashed-dotted line (top left) of the lattice distortion consid- 
ered in the text. The vertical displacement of the carbon 
atoms connected by the hopping integral t± is represented by 
the parameter u. In the top right we schematically represent 
the bandstructure in the vicinity of the Dirac point. 
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plane hopping, t±. From the outset, the atoms lying in 
the A and B sublattices within each layer are not equiv- 
alent, since only one of the sublattices connects to the 
adjacent plane (Fig. 1). This has important experimen- 
tal consequences: one example is the known fact that in 
tunneling experiments one typically detects only one of 
the sublattices of the topmost layer 16 . In the absence of 
a potential difference between the two layers (unbiased 
situation), the bilayer is a zero gap semiconductor, with 
hyperbolic bands touching at the Fermi energy. A change 
in the interlayer hopping will not change this situation, 
and thus the gap can still be tuned through the potential 
difference between layers. However, as shown below, the 
electron-phonon interaction might indeed lead to a stable 
distorted configuration. 

II. THE MODEL 

A. Tight binding description of a biased bilayer 

The electronic Hamiltonian of a biased bilayer consists 
of two contributions, H e = H t b + Hy, where H t b is the 
tight-binding Hamiltonian for the graphene bilayer, and 
Hy reflects the electrostatic bias applied between the 
two graphene planes. The tight-binding Hamiltonian H t b 
contains in itself three terms describing electron itiner- 
ancy among each individual plane and between the two 
planes. In detail we have 

H t b = H t bi + H t b2 + H± , (1) 

with 

H tbl = - *5Z[a^(R)6i ff (R)+a^(R)6i <r (R-ai) 

R,cr 

+ a t la (R)6 1(J (R-a 2 )+H.c.], (2) 
H tb2 = - *5}aUR)MR)+aUR + ai)MR) 

R,cr 

+ 4,(R + a 2 )&2<r(R) + H.c], (3) 
H± = -t ± J2Wl(&)b 2(T (R) + bl(K)a la (K)} , (4) 

R,cr 

and 

Hv= \ 5>L(R)a lCT (R) + &L(R)MR)] 

R,cr 

- \ 5>UR)«2.(R) + &L(n-)MR)] • (5) 

R,(J 

In the above equations ai and a 2 represent the elemen- 
tary translations of the honeycomb lattice. In the pres- 
ence of an electrostatic bias V, the electronic dispersion 
is given by the four branches 

=±^2^ + V 2 + 4* 2 |0 k | 2 ±A k , (6) 



where 

A k = 2^1+4*2(^ + ^2)1^12. (7) 
When V = Eq. (6) simplifies to 

E** = ±\ (± t ± + yJtl+Wlfal*) , (8) 

where 0k is associated with the dispersion of a single 
layer, and is given by 

0k = 1 + e' k ' ai + e' k ' a2 • (9) 
In order to proceed analytically, we approximate |0k| by 

\<t>*\-\aq, (10) 

where we took k = K + q to be close to the Dirac point 
in the honeycomb Brillouin zone, and amounts to using 
the effective mass approximation for bilayer graphene (a 
is the carbon-carbon distance). 

B. Parametrization of Distortion 

We consider a distortion of the perfect lattice struc- 
ture of the bilayer, such that the A and the B atoms, 
connected by the hoping parameter distort along the 
vertical direction by an amount Ui (Fig. 1). In the spirit 
of the SSH model for polyacetylene 14 we assume that, to 
leading order in the strain, the effect of this distortion is 
to change the value of t± according to 

t± = *3_(l + <mi), (11) 

where t \ is the value of the interlayer hoping of the undis- 
torted lattice. For small the in-plane hopping t is 
affected by this distortion only at higher orders in 
and therefore we neglect its variation. This distortion 
will naturally induce an elastic restoring force that we 
parametrize through the term 

N c N c 2 

He^K^uUT^i^ ( 12 ) 

2=1 1=1 

N c denoting the number of unit cells. In the static and 
homogeneous situation the kinetic term gives an average 




FIG. 2: (color online) Phonon modes along the c-axis in 
graphite, using the notation of Ref. 17. The modes A^ u and 
Bi 92 are nearly degenerate at the center of the Brillouin zone, 
with uj ~ 870 cm" 117 ' 18 . 
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null contribution and all U{ acquire the same mean value: 
Ui = u. This leads to a total elastic energy that reads: 



el 



i=l 



N r U 2 



(13) 



The stability analysis of such a distorted phase proceeds 
by minimization of the total electronic and elastic en- 
ergy, given by H e + H e i, with respect to the distortion u. 
We underline that, unlike in the original polyacetylene 
model 14 , the parametrization (11) does not change the 
original periodicity of the lattice, and therefore does not 
require a density commensurability. 
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C. Estimation of parameters 

A precise estimation of the parameters required for 
the computation of the stable distorted configuration is 
not easy nor unique. On the one hand, little is know 
with respect to the structural and elastic properties of 
a graphene bilayer, and thus we will rely on the cor- 
responding knowledge that exists for graphite. On the 
other hand, details like the type of substrate, can signif- 
icantly alter these parameters, as happens, for instance 
with the phonon spectrum that can be sensitive to sub- 
strate and other constraints in the system. 

We will therefore resort to the structural parameters 
(lattice constant and elastic coefficients) known for A- 
B stacked graphite. The carbon-carbon distance is a ~ 
1.42 A, and the graphene unit cell has an area given by 
A c = 3\/3a 2 /2. The equilibrium interlayer distance is 
given by do = co/2 ^ 3.35 A, and corresponds to half the 
unit cell height of A-B stacked graphite 19,20 ' 21 . 

The value of the stiffness, K, can be estimated from 
the phonon spectrum of graphite. In particular the B\ g2 
optical (out-of-plane) phonon mode has a frequency of 
uo ~ 870 cm -1 , which is seen both experimentally 22 , and 
from ab-initio calculations 18 . As a result of the week 
interlayer interaction, this phonon is essentially degener- 
ate with the out-of-plane phonon A2 U present in a single 
layer of graphene. These normal modes are represented 
in Fig. 2. We can assume that K relates to this fre- 
quency through Ka 2 ~ mo; 2 a 2 /4, where m is the car- 
bon atom mass, and a ~ 1.42 A is the carbon-carbon 
distance 23 . As a result we obtain as estimate for the 
stiffness K ~ 8.5eV.A~ 2 . 

With respect to the electron-phonon coupling a, its 
estimation is most straightforward from the knowledge 
of how the interplane hopping varies with distance. 
The interplane hopping, t± corresponds to the tight- 
binding parameter V ppcT in the two center Slater-Koster 
formalism 24,25 . For instance, assuming that 

V ppa (r) ~ Ae~ ar (14) 

one can extract a from interpolation of the hoppings 71 
the in-plane V ppa (a) for graphite. Using the values 17,26 



FIG. 3: (color online) (a) Total energy E* per unit cell as 
a function of the deformation parameter u, and for different 
values of K. The left (right) vertical axis pertains to the solid 
(dashed) curves, (b-d) The equilibrium radius as a function 
of the bias voltage, V. In all panels dashed lines refer to the 
Dirac approximation, whereas full lines have been calculated 
using the full tight-binding dispersion in Eq. (6). Notice that 
in panels (c) and (d) the vertical axis is amplified 10 and 100 
times, respectively. Other parameters used are t° ± = 0.3 eV, 
t 3.0 eV, a 1.5A" 1 . 

71 w 0.4 eV and V ppa (a) w 3.7 eV we obtain a « 1.2 A 1 . 
Alternatively to the formula (14), one could use a more 
refined interpolation formula for V ppcr (r) as discussed in 

Ref. 27. This yields a ~ 1.8 A \ consistent with the 
previous estimate. 

Finally, several recent experiments on the 
bilayer 28,29,30 show that the value of t° ± is essen- 
tially the value expected in graphite, t° ± w 0.3 eV, the 
same applying to the in-plane hopping, t ~ 3 eV. 



III. AB-INITIO CALCULATION OF THE 
ELASTIC CONSTANT 

In addition to the above estimates of the model param- 
eters, we have extracted the compression elastic constant 
from a first principles calculation. Density functional cal- 
culations in graphite and related compounds must be car- 
ried with caution for it is known that different implemen- 
tations of density functional theory can yield noticeably 
different results 18,31 . Having this in mind, we calculated 
the equilibrium distance between graphene planes in the 
bilayer by resorting to two different approximations: the 
Generalized Gradient Approximation (GGA) and the Lo- 
cal Density Approximation (LDA). 

GGA — We sampled the BZ according to the scheme 
proposed by Monkhorst-Pack 32 , with a grid of 12 x 12 x 
4 k-points. Bilayer graphene was modeled in a slab ge- 
ometry by including a vacuum region in a supercell con- 
taining 4 carbon atoms (2 for each graphene sheet). In 
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IV. TOTAL ENERGY AT CONSTANT fi 

Our main objective is to quantify the equilibrium dis- 
tortion that is expected to emerge from the competi- 
tion between elastic and electronic energies in the ground 
state. For illustration purposes we consider first the com- 
putation of the total energy in the (potentially artificial) 
case where the chemical potential is held constant. In 
particular, we assume that the density of carriers in the 
bilayer is such that the chemical potential is located be- 
tween E^Iq and E^Iq: 



V + ^V 2 + 4t 2 ± 



(15) 



FIG. 4: Energy as a function of distance between layers mea- 
sured with respect to the equilibrium position, and calculated 
within the GGA. The positive direction indicates a displace- 
ment towards the other graphene sheet. The zero value energy 
is the energy of the fully relaxed sample. 



Let us start with an unbiased bilayer (V = 0). In this 
case the total electronic energy per unit cell is given by 



The integral is elementary leading to 



(16) 



the normal direction (z-direction), the vacuum separat- 
ing repeating slabs has more than 10 A, and the size of 
the supercell in the z-direction was optimized to make 
sure there was no interaction between repeating slabs. 

LDA — In this case the BZ was sampled within the 
same scheme with a grid of 4 x 4 x 1 k-points, and using 
a supercell comprising 8 carbon atoms (4 for each sheet). 
Adjacent slabs along the z direction were separated by 
more than 30 A, and the size of the supercell along this 
direction was again optimized. 

In either case an increase in the number of sampling 
points did not result in a significant total energy change, 
and the vertical separation quoted above guarantees the 
absence of interaction between adjacent slabs. We used 
dual-space separable pseudopotentials by Hartwigsen, 
Goedecker, and Hutter 33 to describe the ion cores. In 
a first step, all the atoms were fully relaxed to their equi- 
librium positions. Then one of the graphene sheets was 
moved as a whole in the z-direction by very small dis- 
placements, and the total energy of the system was cal- 
culated, without any further relaxation. 

Figure 4 shows the GGA variation of the total energy 
relative to the relaxed sample, as a function of the dis- 
placement from the equilibrium position. Also shown 
is the parabola that was fitted to the calculated values. 
The fitting gave a value of K = 0.615 ± 0.002 eV/A 2 , 
per unit cell. The same calculation within LDA yields 
K = 4.15 ± 0.06 eV/A 2 . The two calculations therefore 
differ by one order of magnitude, signaling the fact that, 
like other systems derived from graphite, density func- 
tional calculations are very sensitive to the details of the 
approximation used. 
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where the momenta qp-, q c and q t are defined as 
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(17) 



(18) 



and A c = 3V3a 2 /2 is the area of the graphene unit cell. 
The total energy E t per unit cell is given by 



(19) 



These two terms compete in such a way that the mini- 
mum energy state is achieved for a finite value of u. 

The dashed lines of the top panels of Fig. 3 represent 
the total energy, E t , as a function of the deformation 
u, using different values of the stiffness parameter, K. 
It is also instructive to investigate to what extent the 
approximation (10) influences the equilibrium deforma- 
tions and, for that, we have performed the calculation 
of the total energies using the full tight-binding disper- 
sions of Eq. (6). The results so obtained are represented 
in the same figure by the solid lines. It is clear from 
Fig. 3(a), that, besides yielding larger slightly larger ab- 
solute values for the energy, the full dispersion increases 
the equilibrium deformation by about 5 to 10%. 

The analytical calculation in the presence of a finite 
bias (V 7^ 0), is also straightforward. The total energy is 
still given by Eq. (16), where E^ ^ is now given by (6). 
The energy integrals are given in the appendix, the final 
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FIG. 5: (color online) Equilibrium deformation, func- 
tion of the electron density per unit cell (n e ) and the bias volt- 
age (V). The parameters used are tj_ = 0.3 eV, t = 3.0 eV, 
K = 0.2 eV.A -2 and a = 1.5 A" 1 . 



instance, through a gate voltage. We define n e as the 
number of electrons per unit cell, with respect to the 
charge-neutral situation in which the valence bands are 
fully occupied. In addition we will be concerned with 
electron doping only. The calculation of the electronic 
energy in this case requires, in general, the consideration 
of three distinct possibilities. Assuming a biased situa- 
tion, and with respect to the notation denned in Fig. 7, 
we can have: 

(i) the Fermi level lying between E\ and E^ in which 
case the Fermi surface consists of a Fermi ring char- 
acterized by two Fermi momenta q^ 1,1 and q^ 1,2 \ 
and the phase space exhibits a central hollow; 

(ii) the Fermi level lying between E2 and £3, where we 
have a more conventional Fermi surface; 



result being 
E e \ = gAc 

N r 27T 



(20) 

The primitives F m ' m (k) are calculated in the appendix, 
with the final result 



Vi 



-^log(2v^R + 2 7 x + 2 m ) J, (21) 

and the remaining parameters are defined in (B5). 

Placing the chemical potential again at the midpoint 
between the two conduction bands at q = 0, 



the corresponding Fermi wavevector is 



(22) 



and we obtain the results shown in the lower panel of 
Fig. 3 for the equilibrium radius. When V varies between 
and 1 eV, the equilibrium radius shows a relative vari- 
ation of ~ 15%. In addition, it can be seen that the dif- 
ference between using the Dirac approximation and the 
full tight-binding dispersion is, in accordance with the 
above, essentially a systematic increase in the equilib- 
rium radius. For this reason, henceforth we will restrict 
the discussion to the results obtained within the Dirac 
approximation. 



(iii) the Fermi level lying above the bottom of the up- 
permost band, in which case we have again two 
Fermi momenta, q^ 1 and \ but the phase space 
is now simply connected. 

The boundaries of these regimes can be easily identified 
through the two threshold densities 



— q 2 , and n e 



^i- (24) 



It follows that the total electronic energy is computed as 
E e gA c 



qdq{E-' + + £"•-) 
27T J q niA q 2tt „ 



(25) 



4 V 



qdqE+'+ 



where the integration limits of the last two terms are 
given by (see also Fig. 7 for notation) 




p0.04- 




V. TOTAL ENERGY AT CONSTANT n e 

We consider now the more relevant case of a bilayer 
with constant carrier density, which can be tuned, for 



FIG. 6: (color online) (a) Selected cuts from Fig. 5 at con- 
stant bias V. For clarity, successive curves have been shifted 
vertically by 0.001 in the order of increasing V. (b) The den- 
sities n% and rie* denned in Eq. (24), plotted as a function of 
V for the same parameters used in Fig. 5. 
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Minimizing the total energy with respect to u yields 
the equilibrium displacements plotted in Fig. 5, for dif- 
ferent electron densities and bias voltages. The typical 
deformations for the parameters quoted in the figure are 
~ 0.11 A, which represents ~ 8% of the carbon-carbon 
distance, a. The variation of u ec[ with n e and V is non- 
monotonic. In particular, one notices that, for constant 
V, the equilibrium deformation tends to saturate beyond 
a given density. This can be appreciated in more detail 
in Fig. 6(a), where we present selected cuts of the same 
surface. The saturation can be understood from the in- 
terplay of two factors: on the one hand, the variation of 
V and n e induce changes in the bandstructure only in a 
region close to the neutrality point; on the other hand, 
for high enough density, the Fermi level will always be 
considerably above the bottom of the uppermost band 
(E3 in Fig. 7). In fact, comparing the values of n* and 
n** presented in Fig. 6(b), one can verify that the first 
sets the scale for the minimum in the curves of u eq _ versus 
n s , therefore defining the shape of the valley in the plot 
of Fig. 5. The value n**, on the other hand, marks the 
onset of saturation. 



VI. DISCUSSION AND CONCLUSIONS 

We have shown that a graphene bilayer with A-B stack- 
ing can be unstable with respect to a Peierls-like distor- 
tion affecting the interplane bonds. This distortion pre- 
serves the bandstructure of the system, in the sense that, 
unlike the original Peierls problem, it does not lead to a 
gap in the unbiased case, nor to its closing in the bi- 



ased situation. In addition, it was found that the general 
effect of the bias voltage is to increase the equilibrium 
deformation. 

By comparing the results obtained with the full tight- 
binding dispersion of bilayer graphene with the effective 
mass approximation [Fig. 3(c-e)], we concluded that the 
former does not introduce significant changes in the equi- 
librium results, and therefore the low energy approxima- 
tion is adequate to study this instability. 

For the values of K used in Fig. 5, the magnitude of 
the deformation corresponds to roughly 10% of the in- 
plane carbon-carbon distance, and is significant. How- 
ever, at this point one can hardly be definite about a 
specific value of the equilibrium deformation on account 
of the uncertainties in the estimation of the parameters 
K and a. The value used for K is close to the compres- 
sive stiffness found with the GGA calculation described 
above. But clearly, had we used the estimate for the 
phonon B\ g2 (or the LDA result) instead, we would have 
obtained much smaller values of i£ e q, as can be inferred 
from Fig. 3(d), although the qualitative features of Fig. 5 
would be preserved. Hence a definitive conclusion as to 
the magnitude of the effect is deferred until the relevant 
parameters in bilayer graphene are experimentally avail- 
able. 

In the consideration of the electronic energy, he have 
accounted only for nearest neighbor in-plane and inter- 
plane hoppings. Additional hopping terms, like next- 
nearest neighbor and other interplane hoppings, should 
not change the qualitative picture presented here. On a 
quantitative level, even the ones that are affected in first 
order in u are expected to contribute only slightly on ac- 
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count of their smaller magnitudes in comparison with t 
and t±. 



Acknowledgments 

V.M.P. is supported by Fundacao para a Ciencia e 
a Tecnologia (FCT) via SFRH/BPD/27182/2006, and 
acknowledges Centro de Fisica do Porto for compu- 
tational support. N.M.R.P and V.M.P acknowledge 
the support of POCI 2010 via PTDC/FIS/64404/2006. 
R.M.R. acknowledges the support of FCT under the 
SeARCH (Services and Advanced Research Comput- 
ing with HTC/HPC clusters) project (contract CONC- 
REEQ/443/2005). 





APPENDIX A: BANDSTRUCTURE 
PARAMETERS 

With respect to the bandstructure depicted in Fig. 7, 
the notable momenta are (vp = 3ta/2): 



Qi 



1 V 4 + 2V 2 t 2 



2v F V t\ + V 2 ' 



V 

Q2 = — , 
Vp v 



while the corresponding energies are 



V 



The energy gap is given by 



A = 2E X 



E 2 = —. 

2 2 ' 

E* = ljV 2 + 4t 2 , 



t±V 



(Ala) 

(Alb) 
(Ale) 

(A2a) 

(A2b) 
(A2c) 

(A3) 



and the midpoint between the upper bands at q = is at 
V+^/V 2 + 4tl 



mid 



(A4) 



to which corresponds the momentum 

q Emid = ^\Jv 2 + 4i? mid + 2^4£; mid (t2 + V*) - t\V* . 

(A5) 



FIG. 7: (color online) Schematic representation of an arbi- 
trary cut of the bandstructure of bilayer graphene close to 
the Dirac point. 



APPENDIX B: ENERGY INTEGRALS 



To compute the total electronic energy, the evaluation 
of the integral 



F ±:± = ±i J kdk\j A + Bk 2 ± 2^D + (Bl) 



is required. The parameters, with respect to the disper- 
sion of the bilayer in Eq. (6), are given as 

A = V 2 + 2t\ , B = Av 2 F , D = t\, 

E = 4v 2 F (tp 2 + V 2 ). (B2) 



The integral is readily computed by changing to the vari- 
able x = VD + Ek 2 , after which it becomes 



+ f3x + jx 2 , 



(B3) 



and is readily available in standard tables. The final 
result is thus 



4v 2 F (t 2 ± + V 2 )\ 3 7 2 7 2 m VH 
"|^ 1o s( 2 V / tR + 2 7^ + 27?2)|, (B4) 
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where 



R = a + fix + 7X 2 , 



A = 4^7 — (3 = 
a = A 



(t\ + V 2 ) 2 ' 
£L> V 4 + d + 3t 2 V 2 



£ ^ + V 2 

t\+Av 2 F (t\ + V 2 )k 2 , 
1 



' ^ + V 2 ' 

/3 = ±2 . (B5) 
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